We need theses packages
library(readxl) # for reading xl files
library(ggplot2) # for figures
library(CaliCo) # calibration package
library(RobustCalibration) # another calibration package
library(DiceDesign) # package for numerical design of experiments
library(DiceKriging) # package for computing Gaussian process Emulator
library(tidyr) # manipulate dataWe load these home-made functions that are the computer models to calibrate
We load some data about ball dropping. The data are used in (Silva et al. 2020) can be found on the github repo.
# the sheet number correspond to the type of ball (see the xl file)
don=read_xlsx("Ball_drops_data.xls",sheet = 5)
names(don)=c("drop","time","Height","Velocity")
don$drop=as.factor(don$drop)We start with a plot
We’ll keep the data from the first drop.
We consider the model (compute code) that provides the height of the
ball as a function of time \(x=t\) and
with two parameters to calibrate \(\theta=(g,h_0)\) the gravity constant and
the initial height. \[(x,\theta)\mapsto
h=f(x,\theta)\] We set the calibration with CaliCo
under model 1 according to (Carmassi et
al. 2019) terminology:
Call:
[1] "model1"
With the function:
function (t, theta)
{
g = theta[1]
h0 = theta[2]
h = -g * 0.5 * t^2 + h0
h[h < 0] = 0
return(h)
}
No surrogate is selected
No discrepancy is added
# we provide some parameters to have a plot
model1 %<% list(theta=c(9.8,46),var=1e-1)
model1$plot(x=don1$time)We can choose the prior distributions on the two parameters in \(\theta\) and on the observational variance (or measurement error):
opt.prior <- list(c(10,1),c(45,5),c(1e-2,1e-0))
type.prior <- c(rep("gaussian",2),"unif")
pr1 <- prior(type.prior,opt.prior)
plot(pr1$Prior1)We can then run the calibration and have a look at the results:
opt.estim = list(Ngibbs=10000,Nmh=10000,thetaInit=c(9.8,45,1e-1),
r=c(0.1,0.1,.1),sig=(diag(c(0.25,1,1e-1)*.3))^2,Nchains=1,burnIn=2000)
mdfit1 <- calibrate(model1,pr1,opt.estim)Call:
With the function:
function (t, theta)
{
g = theta[1]
h0 = theta[2]
h = -g * 0.5 * t^2 + h0
h[h < 0] = 0
return(h)
}
<bytecode: 0x56168fa094d0>
Selected model : model1
Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "38.08%" "24.03%" "28.6%"
Acceptation rate of the Metropolis Hastings algorithm:
[1] "32.17%"
Maximum a posteriori:
[1] 9.2204626 46.7244840 0.1101157
Mean a posteriori:
[1] 9.1313193 46.5504446 0.1039487
By using the MAP for \(\theta\), we can obtain predictions and compute the MSE.
[1] 0.09271534
In order to test the prediction and the coverage rate, we can use cross validation
mdfitCV <- calibrate(model1,pr1,
opt.estim = list(Ngibbs=1000,
Nmh=5000,
thetaInit=c(9.8,45,1e-1),
r=c(0.1,0.1,.1),sig=(diag(c(0.25,1,1e-1)*.3))^2,Nchains=1,burnIn=2000),
opt.valid = list(type.valid="loo",
nCV=50))Call:
With the function:
function (t, theta)
{
g = theta[1]
h0 = theta[2]
h = -g * 0.5 * t^2 + h0
h[h < 0] = 0
return(h)
}
<bytecode: 0x55a30dead190>
Selected model : model1
Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "63.5%" "17.3%" "20.5%"
Acceptation rate of the Metropolis Hastings algorithm:
[1] "59.76%"
Maximum a posteriori:
[1] 9.2636016 46.7784532 0.1069267
Mean a posteriori:
[1] 9.1392385 46.5589027 0.1019072
Cross validation:
Method: loo
Predicted Real Error
1 43.13085 43.13187 0.001017761
2 41.67053 41.89956 0.229031427
3 31.15424 31.54828 0.394039020
4 23.69923 23.83016 0.130929637
5 19.42488 19.44317 0.018283849
6 39.59713 39.91591 0.318785698
RMSE: [1] 0.506053
Cover rate:
[1] "100%"
Predicted Real Error
1 43.13085 43.13187 0.0010177610
2 41.67053 41.89956 0.2290314272
3 31.15424 31.54828 0.3940390202
4 23.69923 23.83016 0.1309296366
5 19.42488 19.44317 0.0182838492
6 39.59713 39.91591 0.3187856982
7 46.55498 46.10805 0.4469310920
8 27.59080 27.89219 0.3013968836
9 16.38248 16.11230 0.2701750090
10 46.52468 46.13628 0.3883948585
11 11.41430 10.83166 0.5826354628
12 32.23738 32.52318 0.2857968186
13 42.26879 42.44439 0.1756022810
14 44.10040 43.89717 0.2032301222
15 34.29827 34.63545 0.3371832491
16 20.85504 20.98674 0.1317023705
17 42.85000 42.93345 0.0834527333
18 17.92037 17.73712 0.1832528099
19 46.39718 46.02302 0.3741640302
20 11.47859 10.83166 0.6469234430
21 44.53033 44.32234 0.2079927645
22 19.42456 19.44317 0.0186116690
23 44.52353 44.32234 0.2011909910
24 33.30193 33.66056 0.3586250895
25 32.23326 32.52318 0.2899216900
26 26.34542 26.51110 0.1656788965
27 42.28382 42.44439 0.1605746881
28 38.82019 39.10350 0.2833092276
29 45.52874 45.28606 0.2426833398
30 42.86560 42.93345 0.0678502795
31 42.85940 42.93345 0.0740591361
32 46.49433 46.05136 0.4429686972
33 27.59683 27.89219 0.2953659406
34 27.59944 27.89219 0.2927532708
35 41.65894 41.89956 0.2406195746
36 38.82299 39.10350 0.2805170876
37 19.43654 19.44317 0.0066320735
38 46.54895 46.10805 0.4408978018
39 36.22515 36.58504 0.3598898127
40 46.56236 46.13628 0.4260795823
41 19.44277 19.44317 0.0003965228
42 13.11581 12.61896 0.4968491463
43 43.63901 43.55704 0.0819716953
44 27.61064 27.89219 0.2815574037
45 32.22386 32.52318 0.2993148857
46 20.91339 20.98674 0.0733543839
47 36.23090 36.58504 0.3541414523
48 44.87945 44.60579 0.2736662923
49 44.89882 44.60579 0.2930376660
50 45.57710 45.28606 0.2910443663
We propose then to add a discrepancy term in the model which corresponds to model 3 according to (Carmassi et al. 2019)’s terminology.
model3 <- model(balldropg,don1$time,don1$Height,"model3",
opt.disc = list(kernel.type="gauss"))
model3 %<% list(theta=c(9.8,46),thetaD=c(1e-1,1),var=1e-3)
plot(model3,don1$time)type.prior <- c(rep("gaussian",2),"unif","unif","unif")
opt.prior <- list(c(10,1),c(40,5),c(1e-2,1e-0),c(0,2),c(1e-2,1e-0))
pr2 <- prior(type.prior,opt.prior)
opt.estim2=list(Ngibbs=1e4,Nmh=1e4,thetaInit=c(10,46,1e-1,0.5,1e-1),
r=c(0.1,0.1),sig=(diag(c(1,1,1e-1,.5,1e-1)*.2))^2,Nchains=1,burnIn=2000)[1] "The computational time might be long to get the output plot"
Call:
With the function:
function (t, theta)
{
g = theta[1]
h0 = theta[2]
h = -g * 0.5 * t^2 + h0
h[h < 0] = 0
return(h)
}
<bytecode: 0x56157fac3c98>
Selected model : model3
Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "65.66%" "45.17%" "51.22%" "53.36%" "54.74%"
Acceptation rate of the Metropolis Hastings algorithm:
[1] "27.8%"
Maximum a posteriori:
[1] 9.5662782 47.3139099 0.1192344 1.3893981 0.1057993
Mean a posteriori:
[1] 9.22878566 46.43685325 0.10479450 0.74960623 0.07172582
By using the MAP for \(\theta\), we can obtain predictions and compute the MSE.
[1] 0.208324
We assume here that the code modeling the ball drop is time-consuming. We make recourse to a Gaussian Process Emulator. We first consider that there is no discrepancy. Hence, we are in the framework of model 2.
# bounds on theta for the GP emulator
binf = c(8,40)
bsup = c(11,50)
# design of experiments and limited number of evaluations to the computer code
ndesign=10
DOE <- maximinSA_LHS(lhsDesign(ndesign,3)$design)$design
DOE <- unscale(DOE,c(min(don1$time),binf),c(max(don1$time),bsup))
Ysim <- balldropg(DOE[1,1],DOE[1,2:3])
for (i in 2:ndesign){Ysim <- c(Ysim,balldropg(DOE[i,1],DOE[i,2:3]))}We provide the design of numerical experiments with the associated evaluations:
model2 <- model(code=NULL,don1$time,don1$Height,"model2",
opt.gp = list(type="matern5_2", DOE=DOE),
opt.sim = list(Ysim=Ysim,DOEsim=DOE))
ParamList <- list(theta=c(9.8,46),var=1e-1)
model2 %<% ParamListand we plot:
opt.estim = list(Ngibbs=1e4,Nmh=1e4,thetaInit=c(9.8,45,1e-1),
r=c(0.1,0.1,.1),sig=(diag(c(0.25,1,1e-1)*.3))^2,Nchains=1,burnIn=2000)
mdfit2 <- calibrate(model2,pr1,opt.estim)Call:
With the function:
NULL
Selected model : model2
Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "48.48%" "25.08%" "28.01%"
Acceptation rate of the Metropolis Hastings algorithm:
[1] "24.13%"
Maximum a posteriori:
[1] 10.15966591 52.44318131 0.02237512
Mean a posteriori:
[1] 9.97885847 47.97409824 0.02184157
We try to make additional evaluations of the model with the purpose to improve the calibration.
model2 <- model(balldropg,don1$time,don1$Height,"model2",
opt.gp = list(type="matern5_2",DOE=NULL),
opt.emul = list(p=2,n.emul=20,binf=binf,bsup=bsup))
newModel2 <- sequentialDesign(model2,pr1,
opt.estim = list(Ngibbs=100,
Nmh=600,
thetaInit=c(9.8,45,1e-1),
r=c(0.1,0.1,.1),sig=(diag(c(0.25,1,1e-1)*.3))^2,
Nchains=1,
burnIn=200),
k=10)We can print the design with added points and make a plot
doe.X
2.2437746 8.328901 46.79794
1.7921889 10.703121 41.68269
1.3998149 8.741023 42.65614
1.0054142 10.681355 45.91142
0.2023622 8.928212 44.82056
2.1777346 10.228825 48.02350
0.7633450 10.035511 40.85691
2.6758861 9.622211 43.49616
0.4874267 8.293524 49.29969
1.2030784 9.237339 47.59909
y.new 0.0000000 9.585993 47.48981
y.new 2.7716651 8.771090 46.12897
y.new 0.6683333 8.913166 46.53728
y.new 1.9033322 8.926643 46.61900
y.new 2.4366653 9.025496 46.42469
y.new 0.0000000 8.939712 46.07890
y.new 2.7716651 9.072577 46.50756
y.new 0.2683333 9.072448 46.62499
y.new 1.3683325 9.012115 46.27818
y.new 1.0349994 8.960061 46.26790
y.new 2.0366655 9.119028 46.61820
And then rerun the calibration with the new design.
model2new = model(code=NULL,don1$time,don1$Height,"model2",
opt.gp = list(type="matern5_2", DOE=newModel2$doe.new),
opt.sim = list(Ysim=as.vector(newModel2$GP.new@y),DOEsim=newModel2$doe.new))
mdfit2seq <- calibrate(model2new,pr1,opt.estim)We are considering model 4. But the calibration is really long…
model4 <- model(code=NULL,don1$time,don1$Height,"model4",
opt.gp = list(type="matern5_2", DOE=NULL),
opt.sim = list(Ysim=Ysim,DOEsim=DOE),
opt.disc = list(kernel.type="gauss"))
optimisation start
------------------
* estimation method : MLE
* optimisation method : BFGS
* analytical gradient : used
* trend model : ~1
* covariance model :
- type : matern5_2
- nugget : NO
- parameters lower bounds : 1e-10 1e-10 1e-10
- parameters upper bounds : 4.975254 5.216081 19.53854
- best initial criterion value(s) : -30.49232
N = 3, M = 5 machine precision = 2.22045e-16
At X0, 0 variables are exactly at the bounds
At iterate 0 f= 30.492 |proj g|= 1.2171
At iterate 1 f = 30.234 |proj g|= 0.10773
At iterate 2 f = 30.197 |proj g|= 0.10669
At iterate 3 f = 30.181 |proj g|= 0.10273
At iterate 4 f = 30.061 |proj g|= 0.49083
At iterate 5 f = 29.969 |proj g|= 0.55312
At iterate 6 f = 29.78 |proj g|= 0.25798
At iterate 7 f = 29.762 |proj g|= 0.072688
At iterate 8 f = 29.76 |proj g|= 0.00743
At iterate 9 f = 29.76 |proj g|= 0.00024236
At iterate 10 f = 29.76 |proj g|= 8.4076e-07
iterations 10
function evaluations 12
segments explored during Cauchy searches 12
BFGS updates skipped 0
active bounds at final generalized Cauchy point 2
norm of the final projected gradient 8.40759e-07
final function value 29.7603
F = 29.7603
final value 29.760286
converged
[1] "The surrogate has been set up, you can now use the function"
Call:
With the function:
NULL
Selected model : model4
Acceptation rate of the Metropolis within Gibbs algorithm:
[1] "71.63%" "58.71%" "65.06%" "67.5%" "68.54%"
Acceptation rate of the Metropolis Hastings algorithm:
[1] "30.92%"
Maximum a posteriori:
[1] 9.3923176 47.7115764 0.1013969 0.6533405 0.1018994
Mean a posteriori:
[1] 9.10329317 46.46869684 0.08681540 0.53280874 0.07886343
[1] "The computational time might be long to get the output plot"
We use the package RobustCalibration that implements the
method proposed in (Gu and Wang 2017)
First, we consider a standard GP for the discrepancy function.
model_gasp=rcalibration(design=don1$time, observations=don1$Height, p_theta=2,simul_type=1,
math_model=balldropg,theta_range=matrix(c(8,11,40,50),2,2,byrow = TRUE)
,S=10000,S_0=2000,discrepancy_type='GaSP')2500 of 10000 posterior samples are drawn
5000 of 10000 posterior samples are drawn
7500 of 10000 posterior samples are drawn
10000 of 10000 posterior samples are drawn
Done with posterior sampling
3759 of 10000 proposed posterior samples of calibration parameters are accepted
5128 of 10000 proposed posterior samples of range and nugget parameters are accepted
m_gasp_pred=predict.rcalibration(model_gasp,don1$time,math_model=balldropg,interval_est=c(0.025,0.975),interval_data = TRUE)25 percent is completed
50 percent is completed
75 percent is completed
100 percent is completed
dfgasp = data.frame(time=don1$time,pred=m_gasp_pred@mean,lowIC=m_gasp_pred@interval[,1],
upIC=m_gasp_pred@interval[,2],obs=don1$Height,
predmathmod=m_gasp_pred@math_model_mean)
dfgasp2 = pivot_longer(dfgasp,cols=2:6,names_to = "type",values_to = "Height")
ggplot()+
geom_line(data = dfgasp2,aes(x=time,y=Height,col=type,linetype=type))+geom_point(data=don1,aes(x=time,y=Height))We also compute the MSEs for the prediction from the mathematical model or from the mathematical model corrected by the discrepancy function.
## MSE if the math model and discrepancy are used for prediction
mean((don1$Height-m_gasp_pred@mean)^2)[1] 0.0007628402
[1] 0.1791588
Second, we consider a scaled GP for the discrepancy function
model_sgasp=rcalibration(design=don1$time, observations=don1$Height, p_theta=2,simul_type=1,
math_model=balldropg,theta_range=matrix(c(8,11,40,50),2,2,byrow = TRUE)
,S=10000,S_0=2000,discrepancy_type='S-GaSP')2500 of 10000 posterior samples are drawn
5000 of 10000 posterior samples are drawn
7500 of 10000 posterior samples are drawn
10000 of 10000 posterior samples are drawn
Done with posterior sampling
2014 of 10000 proposed posterior samples of calibration parameters are accepted
5611 of 10000 proposed posterior samples of range and nugget parameters are accepted
m_sgasp_pred=predict.rcalibration(model_sgasp,don1$time,math_model=balldropg,interval_est=c(0.025,0.975),interval_data = TRUE)25 percent is completed
50 percent is completed
75 percent is completed
100 percent is completed
dfsgasp = data.frame(time=don1$time,pred=m_sgasp_pred@mean,lowIC=m_sgasp_pred@interval[,1],
upIC=m_sgasp_pred@interval[,2],obs=don1$Height,
predmathmod=m_sgasp_pred@math_model_mean)
dfsgasp2 = pivot_longer(dfsgasp,cols=2:6,names_to = "type",values_to = "Height")
ggplot()+
geom_line(data = dfsgasp2,aes(x=time,y=Height,col=type,linetype=type))+geom_point(data=don1,aes(x=time,y=Height))We also compute the MSEs for the prediction from the mathematical model or from the mathematical model corrected by the discrepancy function.
## MSE if the math model and discrepancy are used for prediction
mean((don1$Height-m_sgasp_pred@mean)^2)[1] 0.0008628687
[1] 0.09276917
balldropgfrot that considers an
additional fluid drag force.